#!/usr/bin/python -itt

######### standard library imports ###########################################
import struct
import sys
import time 
from pprint import pprint

######### external library imports ###########################################
import numpy as np
import matplotlib.pyplot as plt

import logging, logging.config, os
scriptpath = os.path.dirname(os.path.realpath(__file__))
logging.config.fileConfig(scriptpath+'/logging.conf')
logger = logging.getLogger('detektorsim')

class Pulse( object ):
    """ A typical G-APD pulse form
    
    Described by:
    f(t) = bsl + n * A0 * heaviside(t-t0)  * (1-exp(-(t-t0)/tau1) * exp((t-t0)/tau2)
    bsl = baseline
    n = number of photons, herre always 1
    A0 = a kind of asymptotical amplitude 
    t0 = the time of start for this discharge
    tau1, tau2, some time constants ... nearly identical I guess.
    
    I leave bsl and A0 away, since these are easy to add from the outside, which leaves me with
    f(t) = heaviside(t-start)  * (1-exp(-(t-start)/tau1) * exp(-(t-start)/tau2)
    
    """
    
    def __init__(self, start=0., tau1=10.*1e-9, tau2=10.*1e-9, xtalk=0.1 ):
        self.start=start
        self.tau1 = tau1
        self.tau2 = tau2
        self.xtalk = xtalk
        
    def __call__(self, t):
        """ Evaluate this function at the times in t
        
        t : *array like*
        """
        
        t = np.array(t)
        if t.ndim != 1:
            raise Exception("t should have ndim 1")
        
        # abreviations
        start = self.start
        tau1 = self.tau1
        tau2 = self.tau2
        
        t1 = t[t < start]
        t2 = (t[t >= start] - start) * -1.
        result = (1. - np.exp(t2/tau1)) * np.exp(t2/tau2)
        result /= 0.25
        
        # Scale the Pulse ... sometimes with 1, sometimes with 2 .. even sometimes with 3 ...
        # and so on... This is my way of adding optical crosstalk
        result *= np.random.geometric(1.-self.xtalk)
        
        return np.concatenate( (np.zeros(t1.shape), result) )

class FasterPulse( object ):
    """ A typical G-APD pulse form
    
    Described by:
    f(t) = bsl + n * A0 * heaviside(t-t0)  * (1-exp(-(t-t0)/tau1) * exp((t-t0)/tau2)
    bsl = baseline
    n = number of photons, herre always 1
    A0 = a kind of asymptotical amplitude 
    t0 = the time of start for this discharge
    tau1, tau2, some time constants ... nearly identical I guess.
    
    I leave bsl and A0 away, since these are easy to add from the outside, which leaves me with
    f(t) = heaviside(t-start)  * (1-exp(-(t-start)/tau1) * exp(-(t-start)/tau2)
    
    """
    
    def __init__(self, thedict):
        
        sipminfo = thedict['settings']['SipmInfo']
        self.tau1 = sipminfo['pulse_tau1']
        self.tau2 = sipminfo['pulse_tau2']
        self.xtalk = sipminfo['crosstalk_probability']
        self.resolution = sipminfo['pulse_max_time_resolution']
        self.sampling_period = sipminfo['pulse_sampling_period']
        
        self.step = int(self.sampling_period/self.resolution)
        
        x = np.arange( 0., 10.*self.tau2, self.resolution)
        x *= -1
        self.values = (1. - np.exp(x/self.tau1)) * np.exp(x/self.tau2)
        
        # With the current tau1 and tau2 settings the half height is reached 
        # after 1ns already. So lets calculate the integral from that time
        # until 15ns later
        integral_start_time = 1e-9
        integral_stop_time = 16e-9
        integral_start_index =  integral_start_time/self.resolution
        integral_stop_index =  integral_stop_time/self.resolution
        
        current_integral = self.values[integral_start_index:integral_stop_index].sum()
        current_integral *= self.resolution
        
        # This integral should be 121 mV ns  
        
        self.values /= current_integral
        self.values *= 121. * 1e-9
        
        # Now the values are in mV 
        
    def __call__(self, starttime, data, weight):

        step = self.step
            
        start = int(starttime/self.resolution)
        s2 = int(starttime/self.sampling_period)
        
        offset = (s2*step-start)%step
        
        values = self.values[offset::step]
        
        if s2 >= 0:
            datalength = len(data[s2:s2+len(values)])
            values = values[:datalength]
            data[s2:s2+len(values)] += values*weight
        else:
            s2 *= -1
            values = values[s2:]
            datalength = len(data[:len(values)])
            values = values[:datalength]
            data[:len(values)] += values*weight
        
            
            


def measureIntegral( thedict, start=0.*1e-9, stop=150*1e-9, step=5e-10):
    """ measure Pulse integral 
    """
    
    pulse_func = thedict['settings']['SipmInfo']['standard_pulse_function']
    
    x = np.arange(start,stop,step)
    y = np.zeros(x.shape)
    pulse_func(0,y,1)
    
    plt.plot(x,y,'.',label="measureIntegral")
    plt.legend()
    plt.grid()
    
    integral = y.sum() * step   # in mV s
    
    return integral
    
    
def calculateAmplitudeWithCrosstalk( thedict, size=1 ):
    
    xtalk = thedict['settings']['SipmInfo']['crosstalk_probability']
    
    amplitude = np.random.geometric(1.-xtalk, size=size)    # amplitude is measured in "Number Of Photons"
    return amplitude


if __name__ == '__main__':
    
    import settings
    thedict = {}
    thedict = settings.put_settings_into_dict( thedict )

    
    plt.ion()
    
    fig2=plt.figure()
    p = thedict['settings']['SipmInfo']['standard_pulse_function']
    x = np.arange(0,150*1e-9,5e-10)
    y = np.zeros(x.shape)
    
    p(25*1e-9,y,1)
    
    # halbe hoehe wird erreicht bei 26*1e-9 s
    # das entspricht x[52] 
    integral_as_jens_and_fabian_use_it = y[52:52+30].sum()
    print "integral_as_jens_and_fabian_use_it", integral_as_jens_and_fabian_use_it
    
    plt.plot(x,y,'.')

    print "measureIntegral"
    print measureIntegral(thedict)
    
